#include <vector>
#include <chrono>
#include <iostream>
#include <immintrin.h>
#include <random>
#include <bitset>
#include<filesystem>
#include<fstream>

void print_double(double value)
{
    unsigned long long int_value;
    memcpy(&int_value, &value, sizeof(value));
    // 打印二进制表示
    // std::cout << "The binary representation of " << value << " is:\n";
    // std::cout << std::bitset<sizeof(int_value) * CHAR_BIT>(int_value) << '\n';

    // 如果你想要分析IEEE 754格式（假设是64位double）
    // 符号位：第1位（从右边开始数，0表示正数，1表示负数）
    // 指数位：接下来的11位
    // 尾数位：剩下的52位
    // std::cout << "Sign bit: " << std::bitset<1>(int_value >> 63) << '\n';
    // std::cout << "Exponent bits: " << std::bitset<11>((int_value >> 52) & 0x7FF) << '\n';
    // std::cout << "Mantissa bits: " << std::bitset<52>(int_value & 0xFFFFFFFFFFFFF) << '\n';
    std::cout << "binary double : "
              << std::bitset<1>(int_value >> 63) << ' '
              << std::bitset<11>((int_value >> 52) & 0x7FF) << ' '
              << std::bitset<52>(int_value & 0xFFFFFFFFFFFFF) << '\n';
}
std::vector<char> input_char;
size_t read_file_to_memory(std::vector<char> &buffer, const std::string &filepath)
{
    // std::vector<char> buffer;
    std::ifstream file(filepath, std::ios::binary | std::ios::ate); // 以二进制模式打开并移动到文件末尾

    if (!file)
    {
        std::cerr << "无法打开文件: " << filepath << std::endl;
        // return buffer; // 返回一个空的 vector
        return 0;
    }

    // 获取文件大小
    std::streamsize size = file.tellg();

    // 移动到文件开头
    file.seekg(0, std::ios::beg);

    // 分配足够的内存
    buffer.resize(size);

    // 读取文件到内存
    if (!file.read(buffer.data(), size))
    {
        std::cerr << "读取文件时出错: " << filepath << std::endl;
        buffer.clear(); // 如果读取失败，则清除 vector
    }

    // 关闭文件
    file.close();

    printf("file byte size=%d long size=%d\n", size, size / 8);
    return size;
    // return buffer;
}
int global_N;
std::vector<double> dvec;
void read_from_file()
{
    FILE *fp_mtx;
    fp_mtx = fopen("E:\\code\\cpp\\other_git_project\\CSR2\\rajat31.mtx", "r");
    int row_num, col_num, nzz_num;
    fscanf(fp_mtx, "%d %d %d", &row_num, &col_num, &nzz_num);
    std::cout << "row_num = " << row_num << "col_num = " << col_num << "nzz_num = " << nzz_num << std::endl;
    dvec.resize(nzz_num);
    global_N = nzz_num;
    //const int N = global_N;
    for (int i = 0; i < nzz_num; ++i)
    {
        int row, col;
        double nzz;
        fscanf(fp_mtx, "%d %d %lf", &row, &col, &nzz);
        dvec[i] = nzz;
    }
    printf("read finish\n");
}
namespace ffc
{
    using U8 = unsigned char;
    using U64 = unsigned long long int;
    using I64 = long long int;
    using F64 = double;
    using U16 = unsigned short;
    const int VECTOR_SIZE = 1 << 10;// 2K*8=16K
    bool is_init = false;

    char ceillog2[256] = {0, 1, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8};

    class information
    {
        // 8+8+2 = 18byte
    public:
        //unsigned short id;
        unsigned long long int stroageEncounterNum{0};
        unsigned long long int EncounterNum_to_EncodeNum{0};
        unsigned char count_1_sum{0};
        unsigned char ceillog2_1_sum{0};
        void test_print()
        {
            //printf("id=%d,stroageEncounterNum=%16llx,EncounterNum_to_EncodeNum=%16llx,count_1_sum=%d\n", id, stroageEncounterNum, EncounterNum_to_EncodeNum, count_1_sum);
        }
    };
    std::vector<information> information_vec(1 << 16); // 占空间较大
    struct state
    {
        // 33byte
        //__m512i m512_decode_offset;//64byte
        U16 isExistRange0_15[16];
        U8 bw; // 存储每个double的bit数  0-64
        //U8 padding[7];
    };

    void init()
    {
        auto start = std::chrono::high_resolution_clock::now();
        for (int i = 1; i < (1 << 16); i++)
        {
            unsigned long long int stroageEncounterNum = 0;
            unsigned long long int EncounterNum_to_EncodeNum = 0;
            unsigned long long int c = 0;
            for (unsigned long long int j = 0; j < 16; ++j)
            {
                unsigned long long int tmp = (i >> j) & 1; // 0 or 1
                stroageEncounterNum |= ((j << (4 * c)) * tmp);
                EncounterNum_to_EncodeNum |= (c << (4 * j));
                c += tmp;
            }
            //information_vec[i - 1].id = i;
            information_vec[i - 1].stroageEncounterNum = stroageEncounterNum;
            information_vec[i - 1].EncounterNum_to_EncodeNum = EncounterNum_to_EncodeNum;
            information_vec[i - 1].count_1_sum = c;
            information_vec[i - 1].ceillog2_1_sum = ceillog2[c - 1]; // c = 1 -> 16
        }
        auto end = std::chrono::high_resolution_clock::now();
        auto duration = std::chrono::duration_cast<std::chrono::nanoseconds>(end - start).count();
        printf("init duration time = %lf ms\n", duration / 1e6); // 0.706ms

        // for(int i=1;i<=(1<<6);++i){
        //     information_vec[i-1].test_print();
        // }
        // information_vec[65534-1].test_print();
        // information_vec[65535-1].test_print();

        is_init = true;
    }
    // fast floating point lossless compress and decompress
    int compress(void *in, int len, void *out)
    {
        // len表示多少个double值
        if (!is_init)
            init();
        U64 *input = (U64 *)in;
        U8 *output = (U8 *)out;
        int length = 0;
        int NUMBER_OF_GROUP = len / VECTOR_SIZE;
        int rest = len - VECTOR_SIZE * NUMBER_OF_GROUP;
        int input_pos = 0, output_pos = 0;
        for (int I = 0; I < NUMBER_OF_GROUP; ++I)
        {
            // int pos = I * VECTOR_SIZE;
            U16 isExistRange0_15[16] = {0};
            for (int i = 0; i < VECTOR_SIZE; ++i)
            {
                U64 num = *(U64 *)(&input[input_pos + i]);
                for (int j = 0; j < 16; ++j)
                {
                    isExistRange0_15[j] |= (1u << ((num >> (60 - j * 4)) & 15));
                }
            }
            state state;
            for (int j = 0; j < 16; ++j)
            {
                state.isExistRange0_15[j] = isExistRange0_15[j];
            }
            unsigned long long int encount_to_encode[16];
            unsigned long long int stroageEncounterNum_4bit[16];
            unsigned long long int decode_offset = 0;
            __m512i m512_encount_to_encode[16];        // 64*8=512bit
            __m512i m512_stroageEncounterNum_4bit[16]; // 32个zmm寄存器
            unsigned char count_1_sum[16];
            unsigned char ceillog2_1_sum[16];
            unsigned char mask[16];
            unsigned char ceillog2_1_prefixsum[16] = {0};
            for (int i = 0; i < 16; ++i)
            {
                encount_to_encode[i] = information_vec[isExistRange0_15[i] - 1].EncounterNum_to_EncodeNum;
                //stroageEncounterNum_4bit[i] = information_vec[isExistRange0_15[i] - 1].stroageEncounterNum;

                m512_encount_to_encode[i] = _mm512_set1_epi64(encount_to_encode[i]);
                //m512_stroageEncounterNum_4bit[i] = _mm512_set1_epi64(stroageEncounterNum_4bit[i]);

                count_1_sum[i] = information_vec[isExistRange0_15[i] - 1].count_1_sum;
                //decode_offset += (count_1_sum[i] == 1) * (((stroageEncounterNum_4bit[i]) & 15ull) << (60 - i * 4));
                ceillog2_1_sum[i] = information_vec[isExistRange0_15[i] - 1].ceillog2_1_sum;
                //mask[i] = ((1 << ceillog2_1_sum[i]) - 1);
                if (i > 0)
                    ceillog2_1_prefixsum[i] = ceillog2_1_sum[i - 1] +  ceillog2_1_prefixsum[i - 1];

                if (0)
                    printf("isExistRange0_15[%d]=%d,encount_to_encode[%d]=%llx,stroageEncounterNum_4bit[%d]=%llx,count_1_sum=%d,ceillog2_1_sum=%d,ceillog2_1_prefixsum=%d\n",
                           i, isExistRange0_15[i] - 1,
                           i, encount_to_encode[i],
                           i, stroageEncounterNum_4bit[i],
                           information_vec[isExistRange0_15[i] - 1].count_1_sum,
                           information_vec[isExistRange0_15[i] - 1].ceillog2_1_sum,
                           ceillog2_1_prefixsum[i]);
            }
            // printf("decode_offset=%llu\n",decode_offset);
            // state.bw = ceillog2_1_prefixsum[15];
            //__m512i m512_decode_offset = _mm512_set1_epi64(decode_offset);
            // state.m512_decode_offset = m512_decode_offset;

            memcpy(&output[output_pos], &state, sizeof(state));
            output_pos += sizeof(state);
            for (int i = 0; i < (VECTOR_SIZE / 8); ++i)
            {
                __m512i m512_input8double = _mm512_loadu_epi64((&input[input_pos + 8 * i]));

                __m512i m512_encode_result = _mm512_set1_epi64(0);
                // m512_encode_result =
                //     ((((m512_encount_to_encode[0] >> (4 * ((m512_input8double >> (60 - 0 * 4)) & 15))) & 15) & mask[0]) << ceillog2_1_prefixsum[0]) +
                //     ((((m512_encount_to_encode[1] >> (4 * ((m512_input8double >> (60 - 1 * 4)) & 15))) & 15) & mask[1]) << ceillog2_1_prefixsum[1]) +
                //     ((((m512_encount_to_encode[2] >> (4 * ((m512_input8double >> (60 - 2 * 4)) & 15))) & 15) & mask[2]) << ceillog2_1_prefixsum[2]) +
                //     ((((m512_encount_to_encode[3] >> (4 * ((m512_input8double >> (60 - 3 * 4)) & 15))) & 15) & mask[3]) << ceillog2_1_prefixsum[3]) +
                //     ((((m512_encount_to_encode[4] >> (4 * ((m512_input8double >> (60 - 4 * 4)) & 15))) & 15) & mask[4]) << ceillog2_1_prefixsum[4]) +
                //     ((((m512_encount_to_encode[5] >> (4 * ((m512_input8double >> (60 - 5 * 4)) & 15))) & 15) & mask[5]) << ceillog2_1_prefixsum[5]) +
                //     ((((m512_encount_to_encode[6] >> (4 * ((m512_input8double >> (60 - 6 * 4)) & 15))) & 15) & mask[6]) << ceillog2_1_prefixsum[6]) +
                //     ((((m512_encount_to_encode[7] >> (4 * ((m512_input8double >> (60 - 7 * 4)) & 15))) & 15) & mask[7]) << ceillog2_1_prefixsum[7]) +
                //     ((((m512_encount_to_encode[8] >> (4 * ((m512_input8double >> (60 - 8 * 4)) & 15))) & 15) & mask[8]) << ceillog2_1_prefixsum[8]) +
                //     ((((m512_encount_to_encode[9] >> (4 * ((m512_input8double >> (60 - 9 * 4)) & 15))) & 15) & mask[9]) << ceillog2_1_prefixsum[9]) +
                //     ((((m512_encount_to_encode[10] >> (4 * ((m512_input8double >> (60 - 10 * 4)) & 15))) & 15) & mask[10]) << ceillog2_1_prefixsum[10]) +
                //     ((((m512_encount_to_encode[11] >> (4 * ((m512_input8double >> (60 - 11 * 4)) & 15))) & 15) & mask[11]) << ceillog2_1_prefixsum[11]) +
                //     ((((m512_encount_to_encode[12] >> (4 * ((m512_input8double >> (60 - 12 * 4)) & 15))) & 15) & mask[12]) << ceillog2_1_prefixsum[12]) +
                //     ((((m512_encount_to_encode[13] >> (4 * ((m512_input8double >> (60 - 13 * 4)) & 15))) & 15) & mask[13]) << ceillog2_1_prefixsum[13]) +
                //     ((((m512_encount_to_encode[14] >> (4 * ((m512_input8double >> (60 - 14 * 4)) & 15))) & 15) & mask[14]) << ceillog2_1_prefixsum[14]) +
                //     ((((m512_encount_to_encode[15] >> (4 * ((m512_input8double >> (60 - 15 * 4)) & 15))) & 15) & mask[15]) << ceillog2_1_prefixsum[15]);
                for (int j = 0; j < 16; ++j)
                {
                    m512_encode_result |= ((((m512_encount_to_encode[j] >> (4 * ((m512_input8double >> (60 - j * 4)) & 15))) & 15)) << ceillog2_1_prefixsum[j]);
                }
                _mm512_storeu_epi64(&output[output_pos + 64 * i], m512_encode_result);
            }
            input_pos += VECTOR_SIZE;
            output_pos += VECTOR_SIZE * 8;

            // U64 sum=1;
            // for(int i=0; i<16; ++i){
            //     sum *= count_1_sum[i];
            // }
            length += ceillog2_1_prefixsum[15] * VECTOR_SIZE + sizeof(state)*8;
            //length += int(ceil(log2(sum))) * VECTOR_SIZE + sizeof(state)*8;
        }

        if (rest)
        {
        }

        return length; // 压缩后的长度
    }

    void decompress(void *in, int len, void *out)
    {
        U8 *input = (U8 *)in;
        U64 *output = (U64 *)out;
        int NUMBER_OF_GROUP = len / VECTOR_SIZE;
        int rest = len - VECTOR_SIZE * NUMBER_OF_GROUP;
        int input_pos = 0;
        int output_pos = 0;
        for (int I = 0; I < NUMBER_OF_GROUP; ++I)
        {
            state state;
            memcpy(&state, &input[input_pos], sizeof(state));
            input_pos += sizeof(state);
            unsigned long long int encount_to_encode[16];
            unsigned long long int stroageEncounterNum_4bit[16];
            unsigned long long int decode_offset = 0;
            __m512i m512_encount_to_encode[16];        // 64*8=512bit
            __m512i m512_stroageEncounterNum_4bit[16]; // 32个zmm寄存器
            unsigned char count_1_sum[16];
            unsigned char ceillog2_1_sum[16];
            unsigned char mask[16];
            unsigned char ceillog2_1_prefixsum[16] = {0};
            U16 *isExistRange0_15 = state.isExistRange0_15;
            for (int i = 0; i < 16; ++i)
            {
                //encount_to_encode[i] = information_vec[isExistRange0_15[i] - 1].EncounterNum_to_EncodeNum;
                stroageEncounterNum_4bit[i] = information_vec[isExistRange0_15[i] - 1].stroageEncounterNum;

                //m512_encount_to_encode[i] = _mm512_set1_epi64(encount_to_encode[i]);
                m512_stroageEncounterNum_4bit[i] = _mm512_set1_epi64(stroageEncounterNum_4bit[i]);

                count_1_sum[i] = information_vec[isExistRange0_15[i] - 1].count_1_sum;
                decode_offset += ((count_1_sum[i] == 1) * (((stroageEncounterNum_4bit[i]) & 15ull) << (60 - i * 4)));
                ceillog2_1_sum[i] = information_vec[isExistRange0_15[i] - 1].ceillog2_1_sum;
                mask[i] = ((1u << ceillog2_1_sum[i]) - 1);

                if (i > 0)
                    ceillog2_1_prefixsum[i] = ceillog2_1_sum[i - 1] +  ceillog2_1_prefixsum[i - 1];

                if (0)
                    printf("isExistRange0_15[%d]=%d,encount_to_encode[%d]=%llx,stroageEncounterNum_4bit[%d]=%llx,count_1_sum=%d,ceillog2_1_sum=%d,ceillog2_1_prefixsum=%d\n",
                           i, isExistRange0_15[i] - 1,
                           i, encount_to_encode[i],
                           i, stroageEncounterNum_4bit[i],
                           information_vec[isExistRange0_15[i] - 1].count_1_sum,
                           information_vec[isExistRange0_15[i] - 1].ceillog2_1_sum,
                           ceillog2_1_prefixsum[i]);
            }
            // printf("decode_offset=%llu\n",decode_offset);
            // state.bw = ceillog2_1_prefixsum[15];
            __m512i m512_decode_offset = _mm512_set1_epi64(decode_offset);

            for (int i = 0; i < (VECTOR_SIZE / 8); i += 1)
            {

                __m512i m512_encode_result = _mm512_loadu_epi64((&input[input_pos + 64 * i])); // 读取8个64 bit值,当做ull处理
                //__m512i m512_decode_result = _mm512_set1_epi64(0);
                __m512i m512_decode_result = m512_decode_offset;

                for (int j = 0; j < 16; ++j)
                {
                    // 运行时循环中不变量，分支预测正确率几乎100%，如何优化为无分支？
                    if (count_1_sum[j] == 1)
                    {
                        /* do nothing */
                    }
                    else if (count_1_sum[j] == 16)[[likely]]
                    {
                        m512_decode_result |= (((m512_encode_result >> ceillog2_1_prefixsum[j]) & 15) << (60 - j * 4));
                    }
                    else
                    {
                        m512_decode_result |= (((m512_stroageEncounterNum_4bit[j] >> (4 * ((m512_encode_result >> ceillog2_1_prefixsum[j]) & mask[j]))) & 15) << (60 - j * 4));
                    }
                    // m512_decode_result |= (( (m512_stroageEncounterNum_4bit[j] >> (4 * ((m512_encode_result >> ceillog2_1_prefixsum[j]) & mask[j]))) & 15) << (60-j*4));
                }

                for (int j = 0; j < 16; ++j)
                {
                    // m512_decode_result |= (( (m512_stroageEncounterNum_4bit[j] >> (4 * ((m512_encode_result >> ceillog2_1_prefixsum[j]) & mask[j]))) & 15) << (60-j*4));
                }
                _mm512_storeu_epi64(&output[output_pos + 8 * i], m512_decode_result);
            }
            input_pos += VECTOR_SIZE * 8;
            output_pos += VECTOR_SIZE;
        }
        if (rest)
        {
        }
    }
};

int main()
{
    read_from_file();
    //read_file_to_memory(input_char,);
    std::random_device rd;
    std::mt19937_64 gen(rd());

    const int N = 1 << 27;

    std::vector<double> input_vector(N);
    std::vector<unsigned long long int> compress_vector(N * 2);
    std::vector<double> decompress_vector(N);

    for (int i = 0; i < N; ++i)
    {
        //input_vector[i] = gen();
    }

    //compress
    {
        auto start1 = std::chrono::system_clock::now();
        int len=ffc::compress(dvec.data(), global_N, compress_vector.data());//N为double值数量,非字节数，实际字节数为N*8
        auto end1 = std::chrono::system_clock::now();
        auto duration1 = std::chrono::duration_cast<std::chrono::nanoseconds>(end1 - start1).count();
        printf("compress time : %lf ms\n", duration1 / 1e6);
        printf("len=%d\n", len);
        printf("ratio = %lf\n",(double)len/(global_N*64));
    }

    //decompress
    {
        auto start = std::chrono::system_clock::now();
        ffc::decompress(compress_vector.data(), global_N, decompress_vector.data());
        auto end = std::chrono::system_clock::now();
        auto duration = std::chrono::duration_cast<std::chrono::nanoseconds>(end - start).count();
        printf("decompress time : %lf ms\n", duration / 1e6);
    }

    int cmp_result = 1;
    for (int i = 0; i < N; ++i)
    {
        // if(input_vector[i] != decompress_vector[i])
        // {
        //     printf("Error at index %d,input_vector[%d]=%lf,decompress_vector[%d]=%lf\n",i,i,input_vector[i],i,decompress_vector[i]);
        //     print_double(input_vector[i]);
        //     print_double(decompress_vector[i]);
        //     return -1;
        // }
        cmp_result ^= (input_vector[i] == decompress_vector[i]);
    }
    if (cmp_result)
    {
        printf("Test passed\n");
    }
    else
    {
        printf("Test failed\n");
    }

    return 0;
}